
#ifndef __CLUSTERING_H_
#define __CLUSTERING_H_

#include "mesh.h"
#include "uskeleton.h"
#include <vector>
#include <queue>

class Clustering {

public:
   struct global_stack_element {
      
      global_stack_element() :
         bone(0), vertex(0), e_init(0), e_term(0), dist_to_skel(0), 
         dist_to_prim_skin(0), component(0) {}
      
      global_stack_element(int b, int v, int ei, int et, Number ds, Number dp, int c) :
         bone(b), vertex(v), e_init(ei), e_term(et), dist_to_skel(ds), 
         dist_to_prim_skin(dp), component(c) {}
      
      int       bone;
      int       vertex;
      int       e_init;
      int       e_term;
      Number    dist_to_skel;  //< distance to skeleton
      Number    dist_to_prim_skin;  //< distance to primary skin
      int       component;
   };

   struct stack_2_distances {
      
      stack_2_distances() :
         vertex(0), e_init(0), e_term(0), dist_to_xv(0), order_of_stack(0) {}
      
      stack_2_distances(const int v, const int ei, const int et, const Number dxv, const Number o) :
         vertex(v), e_init(ei), e_term(et), dist_to_xv(dxv), order_of_stack(o) {}
         
      int        vertex;
      int        e_init;
      int        e_term;
      Number      dist_to_xv;  // distance to the last intersection found
      Number      order_of_stack;
   };

   struct stack_element {
      
      stack_element() :
         vertex(0), e_init(0), e_term(0), dist(0) {}
      
      stack_element(const int v, const int ei, const int et, const Number d) :
         vertex(v), e_init(ei), e_term(et), dist(d) {}
         
      int        vertex;
      int        e_init;
      int        e_term;
      Number      dist;  //< distance to connected component (key)
   };

   struct compare_stackelem {
      bool operator() (const stack_element &a, const stack_element &b) const {
         return a.dist > b.dist;
      }
   };
      
   struct compare_globalstackelem {
      bool operator() (const global_stack_element &a, const global_stack_element &b) const {
         return a.dist_to_skel > b.dist_to_skel;
      }
   };
      
   struct compare_stack2distances {
      bool operator() (const stack_2_distances &a, const stack_2_distances &b) const {
         return a.order_of_stack > b.order_of_stack;
      }
   };

   struct edge_pair {
      
      edge_pair() { e_ini=-1; e_term=-1; };
      int e_ini;
      int e_term;
   };

protected:
   Number length(int id_he) const {
      
      Halfedge_iterator edge = model->halfedges_begin();
      const Point3 &pA = edge[id_he].vertex()->point();
      const Point3 &pB = edge[id_he].opposite()->vertex()->point();
      const Vector3 v = pB - pA;      
      return sqrt(v*v);
   }

   bool interTriSeg(const Vector3 &u, const Vector3 &v, const Vector3 &dir, const Vector3 &w0) {
      
      Vector3  n;             // triangle vectors
      Vector3  w;             // ray vectors
      Vector3  intersection;  // intersection of ray and plane
      Number   r, a, b;       // params to calc ray-plane intersect
      Number   SMALL_NUM = 1E-8;
      
      // get triangle edge vectors and plane normal
      n = cross_product(u,v);             // cross product
      if (n * n == 0)                    // triangle is degenerate
          return false;                 // another triangle will contain the intersection if
                                          //there is one.
      
      a = n*w0;
      b = n*dir;
      if (fabs(b) < SMALL_NUM) {     // ray is parallel to triangle plane
          return false;             // ray disjoint from plane
      }

      // get intersect point of ray with triangle plane
      r = a/b;
      if (r < 0.0 || r > 1.0)                   // segment not intersect triangle plane
          return false;

      intersection = w0 - r*dir;           // intersect point of ray and plane
      // is I inside T?
      Number t1, t2;
      float Det = u[0]*v[1] - u[1]*v[0]; 
      if (fabs(Det) > SMALL_NUM) {
         t1 = (intersection[0]*v[1] - intersection[1]*v[0])/Det;
         if(t1 < 0 || t1 > 1) return false;
         t2 = (intersection[1]*u[0] - intersection[0]*u[1])/Det;
         if(t2 < 0 || t2 > 1 || (t1+t2)>1 ) return false;
         return true;
      }
      Det = u[0]*v[2] - u[2]*v[0];
      if (fabs(Det) > SMALL_NUM) {
         t1 = (intersection[0]*v[2] - intersection[2]*v[0])/Det;
         if(t1 < 0 || t1 > 1) return false;
         t2 = (intersection[2]*u[0] - intersection[0]*u[2])/Det;
         if(t2 < 0 || t2 > 1  || (t1+t2)>1 ) return false;
         return true;
      }
      Det = u[1]*v[2] - u[2]*v[1];
      t1 = (intersection[1]*v[2] - intersection[2]*v[1])/Det;
      if(t1 < 0 || t1 > 1) return false;
      t2 = (intersection[2]*u[1] - intersection[1]*u[2])/Det;
      if(t2 < 0 || t2 > 1  || (t1+t2)>1 ) return false;
      return true;
   }
   
protected:
   Mesh        *model;     //the model represented by a mesh and a skeleton
   USkeleton   *skeleton;  //
   
public:
   int    lista_size;  // size of the list used to find a component by horizontal search
};

#endif
